

*=============================================
* Figure B.5 Predicted BEV and ICEV ownership
*=============================================

use    "${dataout}MainDataset" , clear

keep if year>=2015 & year<=2017
keep if couple==1

/*******************************************************************/
/**********************  BEV outcome *******************************/
/*******************************************************************/

preserve
	** Select variables to loop over: 
	local outcomelist  BEV_fam_yes  
	local treatvar1list toll_fam_mean_KPI  
	local treatvar2list ptl_fam_km_mean        

	** Regressions - loop over: outcome, treatment, year
	foreach outcome in `outcomelist'{
		foreach treatvar1 in `treatvar1list'{
			foreach treatvar2 in `treatvar2list'{

				reghdfe `outcome' `treatvar1'  `treatvar2'  $age    $distance  $time  $publictime  $publicquality  , ///
				absorb($FE   $household  $income  $employment  $education  i.year  , savefe ) vce(cluster $clustervar) 
				summarize `outcome' if e(sample)==1
				estadd scalar MeanDep = r(mean), replace
				summarize `treatvar1' if e(sample)==1
				estadd scalar Treat1 = r(mean), replace
				summarize `treatvar2' if e(sample)==1
				estadd scalar Treat2 = r(mean), replace
				estadd local year="2015-17" , replace
				estadd local grkFE="\checkmark", replace
				estadd local grkbFE="\checkmark", replace
				estadd local HHcon="\checkmark", replace
				estadd local HHWork="\checkmark", replace
				estadd local PublicTime="\checkmark", replace
			}
		}
	}


	** Mark main sample
	capt drop MainSample
	gen MainSample=0
	replace MainSample=1 if e(sample)==1

	// NB NB
	keep if MainSample==1 

	display "Mean dep: `e(MeanDep)'"

	gen beta_toll=_b[toll_fam_mean_KPI]
	gen beta_ptl = _b[ptl_fam_km_mean]

	predict yhat_xb, xb
	predict yhat_xbd , xbd   /*xb + d*/
	predict yhat_xbd_se , stdp  

	gen yhat_xbd_toll_0 =yhat_xbd - (beta_toll * toll_fam_mean)
	gen yhat_xbd_toll_0_ptl_0 =yhat_xbd - (beta_toll * toll_fam_mean) - (beta_ptl * ptl_fam_km_mean)

	* Making figure
	graph bar ///
		(mean) yhat_xbd ///
		(mean) yhat_xbd_toll_0 ///
		(mean) yhat_xbd_toll_0_ptl_0  ///
	if MainSample==1 ///
	, graphregion(color(white)) /// 
	legend(cols(1) label(1 "Actual road toll and bus lane") ///
	label(2 "0 road toll and actual bus lane") ///
	label(3 "0 road toll and 0 bus lane")) ytitle("BEV (yes=1)") ///
	blabel(bar, position(bar) format(%9.4f)tstyle(body) color(black))  ///
	bar(1,color(black) fcolor(navy)) ///
	bar(2,color(black) fcolor(navy*0.60)) ///
	bar(3,color(black) fcolor(navy*0.30))  scale(1.3)

	graph export     "${figures}FigureB5a.png" , replace
	graph export     "${figures}FigureB5a.pdf" , replace
	graph save       "${figures}FigureB5a.gph" , replace
restore 

/*******************************************************************/
/**********************  ICE outcome *******************************/
/*******************************************************************/

preserve
	** Select variables to loop over: 
	local outcomelist  ICE_fam_count   
	local treatvar1list toll_fam_mean_KPI  
	local treatvar2list ptl_fam_km_mean        

	** Regressions - loop over: outcome, treatment, year
	foreach outcome in `outcomelist'{
		foreach treatvar1 in `treatvar1list'{
			foreach treatvar2 in `treatvar2list'{
				reghdfe `outcome' `treatvar1'  `treatvar2'  $age    $distance  $time  $publictime  $publicquality  , ///
				absorb($FE   $household  $income  $employment  $education  i.year  , savefe ) vce(cluster $clustervar) 
				summarize `outcome' if e(sample)==1
				estadd scalar MeanDep = r(mean), replace
				summarize `treatvar1' if e(sample)==1
				estadd scalar Treat1 = r(mean), replace
				summarize `treatvar2' if e(sample)==1
				estadd scalar Treat2 = r(mean), replace
				estadd local year="2015-17" , replace
				estadd local grkFE="\checkmark", replace
				estadd local grkbFE="\checkmark", replace
				estadd local HHcon="\checkmark", replace
				estadd local HHWork="\checkmark", replace
				estadd local PublicTime="\checkmark", replace
			}
		}
	}

	** Mark main sample
	capt drop MainSample
	gen MainSample=0
	replace MainSample=1 if e(sample)==1

	* Dropping others
	keep if MainSample==1 

	display "Mean dep: `e(MeanDep)'"

	gen beta_toll=_b[toll_fam_mean_KPI]
	gen beta_ptl = _b[ptl_fam_km_mean]

	predict ICEhat_xbd , xbd   /*xb + d_absorbvars*/
	predict ICEhat_xbd_se , stdp

	gen ICEhat_xbd_toll_0 =ICEhat_xbd - (beta_toll * toll_fam_mean_KPI )
	gen ICEhat_xbd_toll_0_ptl_0 =ICEhat_xbd - (beta_toll * toll_fam_mean_KPI ) - (beta_ptl  * ptl_fam_km_mean )

	sum ICEhat_xbd ICEhat_xbd_toll_0 ICEhat_xbd_toll_0_ptl_0

	* Making bar chart
	graph bar ///
		(mean) ICEhat_xbd ///
		(mean) ICEhat_xbd_toll_0 ///
		(mean) ICEhat_xbd_toll_0_ptl_0  ///
	if MainSample==1 ///
	, graphregion(color(white)) ///
	legend(cols(1) label(1 "Actual road toll and bus lane") ///
	label(2 "0 road toll and actual bus lane") ///
	label(3 "0 road toll and 0 bus lane")) ytitle("Number of ICEVs") ///
	blabel(bar, position(bar) format(%9.4f)tstyle(body) color(black))  ///
	bar(1,color(black) fcolor(navy)) ///
	bar(2,color(black) fcolor(navy*0.60)) ///
	bar(3,color(black) fcolor(navy*0.30)) scale(1.3)

	graph export     "${figures}FigureB5b.png" , replace
	graph export     "${figures}FigureB5b.pdf" , replace
	graph save       "${figures}FigureB5b.gph" , replace
restore
